Generalidades

En esta entrega se hace una verificación complementaria en función de:

• Producción de sedimentos: Se muestra la producción anual de sedimentos para varias tipos de uso del suelo y pendientes.
• Cantidad de agua: Se incluyen estaciones de caudal que fueron descartadas en el proceso de calibración/validación.

Sedimentos

La siguiente figura muestra un resumen de la producción anual de sedimentos ordenado según el uso de la tierra. Los usos EUCA y MONT tienen una producción de sedimentos bastante inferior AGRP, AGRC, LECH y GRAS. Adicionalmente AGRP, LECH, y GRAS tienen un comportamiento bastante similar. AGRC está ligeramente por encima y su variabilidad practicamente duplica la variabilidad del grupo AGRP, LECH y GRAS.

SYLD <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/SYLD_promedio_10anos.RDS")
hru_info <- readRDS("~/R_projects/SWAT_calibration_ETmodis/shapefiles_RDS/hru_info.RDS")

SYLD_GRAS = unlist(SYLD[(hru_info$LANDUSE=="GRAS")])
SYLD_EUCA = unlist(SYLD[(hru_info$LANDUSE=="EUCA")])
SYLD_MONT = unlist(SYLD[(hru_info$LANDUSE=="MONT")])
SYLD_AGRC = unlist(SYLD[(hru_info$LANDUSE=="AGRC")])
SYLD_AGRP = unlist(SYLD[(hru_info$LANDUSE=="AGRP")])
SYLD_LECH = unlist(SYLD[(hru_info$LANDUSE=="LECH")])

library(plotly)
fig <- plot_ly(y = ~SYLD_AGRC, type = "box", name="AGRC")
fig <- fig %>% add_trace(y = ~SYLD_AGRP, name="AGRP")
fig <- fig %>% add_trace(y = ~SYLD_LECH, name="LECH")
fig <- fig %>% add_trace(y = ~SYLD_GRAS, name="GRAS")
fig <- fig %>% add_trace(y = ~SYLD_EUCA, name="EUCA")
fig <- fig %>% add_trace(y = ~SYLD_MONT, name="MONT")
fig <- fig %>% layout(yaxis = list(title="Producción anual de sedimentos [ton/ha]"), 
                      showlegend = FALSE,
    autosize=FALSE,
    width=500,
    height=350)
fig

A continuación se muestra la producción de sedimentos en funcion de las 6 clases de pendientes modeladas en SWAT:

Sband = paste(unique(hru_info$SLOPE_BAND))
paste0("S", 1:6,"(",Sband,")")

## [1] "S1(0-3.0)"     "S2(3.0-6.0)"   "S3(6.0-9.0)"   "S4(9.0-12.0)" 
## [5] "S5(12.0-15.0)" "S6(15.0-9999)"

slope_plot = function(Luse){
c1 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[1]
S1 = unlist(SYLD[c1])
c2 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[2]
S2 = unlist(SYLD[c2])
c3 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[3]
S3 = unlist(SYLD[c3])
c4 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[4]
S4 = unlist(SYLD[c4])
c5 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[5]
S5 = unlist(SYLD[c5])
c6 = hru_info$LANDUSE==Luse & hru_info$SLOPE_BAND==Sband[6]
S6 = unlist(SYLD[c6])

if(is.null(S1)){S1=0}
if(is.null(S2)){S2=0}
if(is.null(S3)){S3=0}
if(is.null(S4)){S4=0}
if(is.null(S5)){S5=0}
if(is.null(S6)){S6=0}

fig <- plot_ly(y = ~S1, type = "box", name=paste0("s1 (",table(c1)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S2, name=paste0("s2 (",table(c2)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S3, name=paste0("s3 (",table(c3)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S4, name=paste0("s4 (",table(c4)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S5, name=paste0("s5 (",table(c5)["TRUE"],")" ))
fig <- fig %>% add_trace(y = ~S6, name=paste0("s6 (",table(c6)["TRUE"],")" ))
fig <- fig %>% layout(title=paste(Luse),
                      yaxis = list(title="Producción anual de sedimentos [ton/ha]",
                                   range = c(0, 70)), 
                      showlegend = FALSE,
                      autosize=FALSE,
                      width=350,
                      height=300)
fig
}


leafsync::sync(slope_plot("AGRC"),slope_plot("AGRP"),
             slope_plot("LECH"),slope_plot("GRAS"))

Cantidad de agua

La verificación de cantidad de agua se hace a escala mensual en el período 2006-2015 (2005 de calentamiento). Esta verificacion muetra a grandes trazos que en los primeros años de simulacion (2006-2010) predomina la subestimación. Luego a partir del año 2010 se sobreestima el caudal.

A continuación leemos los resultados del modelo ejecutado en ese período así como los valores observados. Nótese que la simulación de caudal tiene de 100 ejecuciones del modelo, los estadísticos se dan en funcion de la mediana de las 100 simulaciones.

flow_m23 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/flow_m23.RDS")
flow_m30 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/flow_m30.RDS")
flow_m36 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/flow_m36.RDS")
flow_m38 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/flow_m38.RDS")
flow_m41 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/optim20200401_22/flow_m41.RDS")

flow_m23 = lapply(flow_m23[,-1],function(x){convertFlow(zoo(x, flow_m23$date),
                                                        from="cumecs", to="mm", area.km2 = 687,
                                                        timestep.default = "months")}) 
flow_m30 = lapply(flow_m30[,-1],function(x){convertFlow(zoo(x, flow_m30$date),
                                                        from="cumecs", to="mm", area.km2 = 1092,
                                                        timestep.default = "months")}) 
flow_m36 = lapply(flow_m36[,-1],function(x){convertFlow(zoo(x, flow_m36$date),
                                                        from="cumecs", to="mm", area.km2 = 2744,
                                                        timestep.default = "months")})
flow_m38 = lapply(flow_m38[,-1],function(x){convertFlow(zoo(x, flow_m38$date),
                                                        from="cumecs", to="mm", area.km2 = 3159,
                                                        timestep.default = "months")}) 
flow_m41 = lapply(flow_m41[,-1],function(x){convertFlow(zoo(x, flow_m41$date),
                                                        from="cumecs", to="mm", area.km2 = 4896,
                                                        timestep.default = "months")})

Qsim = cbind(Qsim_23 = apply(as.data.frame(flow_m23),1,median),
             Qsim_30 = apply(as.data.frame(flow_m30),1,median),
             Qsim_36 = apply(as.data.frame(flow_m36),1,median),
             Qsim_38 = apply(as.data.frame(flow_m38),1,median),
             Qsim_41 = apply(as.data.frame(flow_m41),1,median))


Qobs <- readRDS("~/R_projects/SWAT_calibration_ETmodis/data/runoff_mensual_2006_2015.RDS")

qo = cbind(Qobs, Qsim)

La siguiente tabla muestra los estadísticos KGE, NSE y Pbias a escala mensual en el período 2006-2015.

library(hydroGOF)

kge_stat = round(c(KGE(qo$Qsim_23,qo$Qobs_23),
             KGE(qo$Qsim_30,qo$Qobs_30),
             KGE(qo$Qsim_36,qo$Qobs_36),
             KGE(qo$Qsim_38,qo$Qobs_38),
             KGE(qo$Qsim_41,qo$Qobs_41)),2)

nse_stat = round(c(NSE(qo$Qsim_23,qo$Qobs_23),
             NSE(qo$Qsim_30,qo$Qobs_30),
             NSE(qo$Qsim_36,qo$Qobs_36),
             NSE(qo$Qsim_38,qo$Qobs_38),
             NSE(qo$Qsim_41,qo$Qobs_41)),2)

pbias_stat = round(c(pbias(qo$Qsim_23,qo$Qobs_23),
             pbias(qo$Qsim_30,qo$Qobs_30),
             pbias(qo$Qsim_36,qo$Qobs_36),
             pbias(qo$Qsim_38,qo$Qobs_38),
             pbias(qo$Qsim_41,qo$Qobs_41)),1)

stat = rbind(KGE=kge_stat,NSE=nse_stat, PBIAS=pbias_stat)

kable(stat, col.names = c("WL_23", "WL_30", "WL_36", "WL_38", "WL_41")) %>%
  kable_styling(c("striped", "bordered"))

res = cbind(r23 = (qo$Qsim_23 - qo$Qobs_23),
            r30 = (qo$Qsim_30 - qo$Qobs_30),
            r36 = (qo$Qsim_36 - qo$Qobs_36),
            r38 = (qo$Qsim_38 - qo$Qobs_38),
            r41 = (qo$Qsim_41 - qo$Qobs_41))

Res_plot = function(x){
  
  cc = substr(index(res),1,3) %in% "ene"
  nn = index(res)
  nn[!cc] = NA 
  nn = substr(nn,6,10)
  
  barplot(res[,x], las=1, col=ifelse(res[,x]>=0,4,2),
          main = x, ylab="Qsim-Qobs (mm)", names=nn, ylim = c(-100,100))
} 

par(mfrow=c(5,1), mar=c(4,4,1,1))
Res_plot("r23")
Res_plot("r30")
Res_plot("r36")
Res_plot("r38")
Res_plot("r41")